Analyzing packets captured on 900MHz ISM

I captured a few thousand(!) packets in 12 seconds(!!) of snooping the entire 900MHz ISM band. These were extracted in code parallel to the ISM Snooping example, with the final bitstream state dumped out in the JSON format at the end. We'll now pull that in and look at the bitstreams.

In [154]: 
%pylab inline

pylab.rcParams["figure.figsize"] = 14,4

import scipy
import scipy.signal

import gc

import os
import json

Populating the interactive namespace from numpy and matplotlib


WARNING: pylab import has clobbered these variables: ['f']
`%matplotlib` prevents importing * from pylab and numpy
In [155]: 
FILENAME = "capture_state_1421785415.json"

f = file(FILENAME)
capture_state = json.load(f)
f.close()

capture_state.keys()

ERROR! Session/line number was not unique in database. History logging moved to new session 168
Out[155]:
[u'merged_hits',
 u'Fs',
 u'packets',
 u'filename',
 u'sample_size',
 u'Fc',
 u'raw_hits']
In [156]: 
[ capture_state[x] for x in ["filename", "Fc", "Fs"] ]
Out[156]:
[u'/home/jbm/grc/all_of_900_ism.sdr', 915000000.0, 26000000.0]
In [157]: 
capture_state['packets'][0].keys()
Out[157]:
[u'bitstream', u'hit_i', u'baudrate']
In [158]: 
p = capture_state["packets"][1]
plot(p['bitstream'])
ylim(-0.1, 1.1)
title("Hit %d, baudrate is %d" % (p["hit_i"], p["baudrate"]))
Out[158]:
<matplotlib.text.Text at 0x7fb467ff3ad0>
In [159]: 
packets = [ p for p in capture_state["packets"] if len(p['bitstream']) > 0 ]
hits = capture_state["merged_hits"]
In [160]: 
len(packets)
Out[160]:
2339
In [161]: 
hist([p['baudrate'] for p in packets ], bins=100); None
In [162]: 
hist([np.log10(len(p['bitstream'])) for p in packets], bins=100); None
In [163]: 
packets_peakish = [ p for p in packets if len(p['bitstream']) > 100 and len(p['bitstream']) < 150 ]
hist([len(p['bitstream']) for p in packets_peakish], bins=range(100,150)); None
In [164]: 
packets_peakish = [ p for p in packets if len(p['bitstream']) > 115 and len(p['bitstream']) < 135 ]
hist([len(p['bitstream']) for p in packets_peakish], bins=range(100,150)); None
In [165]: 
hist([(p['baudrate']) for p in packets_peakish], bins=100); None
In [166]: 
packets_peakish_slow = [ p for p in packets_peakish if p['baudrate'] < 20000]
hist([len(p['bitstream']) for p in packets_peakish_slow], bins=range(100,150)); None
In [167]: 
to_show = []

for p in packets_peakish_slow:
    to_add = p['bitstream'] + [0] * (125 - len(p['bitstream']))
    to_show.append(to_add)
    
to_show = np.array(to_show)
to_show.shape
Out[167]:
(224, 125)
In [168]: 
imshow(to_show)
Out[168]:
<matplotlib.image.AxesImage at 0x7fb467e4f9d0>
In [169]: 
def bitlen_rle(b):
    "Returns a RLE of the central binary state lengths, with no concept of parity. Ignores ends."
    edges = (b[:-1] != b[1:]).nonzero()[0]
    # XXX This neglects the first and last segments. We don't really care, but.
    run_lengths = np.diff(edges)
    return run_lengths

def find_preamble_end(bitstream):
    bp = bitstream[:-1] != bitstream[1:]
    bpr = bitlen_rle(bp)
    bitpos = 0
    
    for j in bpr:
        if j >= 8:
            return bitpos+j
        bitpos += j
        
    return None

def find_preamble_perfectly(bitstream):
    starting_point = find_preamble_end(bitstream)
    if None == starting_point:
        return None
    starting_point -= 5
    if starting_point < 0:
        starting_point = 0
        
    last_state = bitstream[starting_point]
    
    c = starting_point+1
    while c < len(bitstream):
        if bitstream[c] == last_state:
            return c
        last_state = bitstream[c]
        c += 1        
    return c
In [170]: 
to_show = []

for p in packets_peakish_slow:
    preamble_header = find_preamble_perfectly(np.array(p['bitstream']))
    if None == preamble_header:
        continue # preamble_header = 0
    to_add = p['bitstream'][(preamble_header):] + [0] * (125 - len(p['bitstream']) + preamble_header)
    to_show.append(to_add)

to_show.sort()
to_show = np.array(to_show)
to_show.shape
Out[170]:
(218, 125)
In [171]: 
imshow(to_show, cmap=cm.gray_r, interpolation="none")
Out[171]:
<matplotlib.image.AxesImage at 0x7fb45f5e5bd0>
In [172]: 
to_show[100,0:25]
Out[172]:
array([1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1,
       1, 1])
In [173]: 
preamble_header = [1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1]
ph = "".join(map(str, preamble_header))

def find_actual_preamble(bitstream):
    bs = "".join(map(str, bitstream))
    
    return bs.find(ph)
In [174]: 
to_show = []

for p in packets_peakish_slow:
    preamble_header = find_actual_preamble(np.array(p['bitstream']))
    if -1 == preamble_header:
        continue # preamble_header = 0
    to_add = p['bitstream'][(preamble_header):] + [0] * (125 - len(p['bitstream']) + preamble_header)
    to_show.append(to_add)

to_show.sort()
to_show = np.array(to_show)
to_show.shape
Out[174]:
(212, 125)
In [175]: 
imshow(to_show[:,30:60].T, cmap=cm.gray_r, interpolation="none")
Out[175]:
<matplotlib.image.AxesImage at 0x7fb45e1bfc10>
In [176]: 
real_data = to_show[:, 32:70]
imshow(real_data.T, interpolation="none", cmap=cm.gray_r)
Out[176]:
<matplotlib.image.AxesImage at 0x7fb45f32e890>

Pretty. How does the bit distribution look per-line?

In [177]: 
plot(np.sum(real_data, axis=1))
Out[177]:
[<matplotlib.lines.Line2D at 0x7fb45f4b6250>]

So there's a variable number of bits set per-entry at this level. How does each bit do across frames/packets?

In [195]: 
real_data = to_show[:, 32:70]

plot(np.sum(real_data, axis=0)/float(real_data.shape[0]))
axhline(0.5, color="red", alpha=0.5)
axvline(3, color="green", alpha=0.7)
axvline(25, color="green", alpha=0.7)

real_data = to_show[:, 32+3:32+25]

Hey now, that looks like our actual packet body. Nice encoding work, people who designed this. Focusing down onto that data range.

In [231]: 
plot(np.sum(real_data, axis=0)/float(real_data.shape[0]))
axhline(0.5, color="red", alpha=0.5)
Out[231]:
<matplotlib.lines.Line2D at 0x7fb45de266d0>
In [232]: 
imshow(real_data.T, interpolation="none", cmap=cm.gray_r)

i = 0
while False: # i < real_data.shape[1]:
    axhline(i, color="red")
    i += 8
In [233]: 
plot(np.sum(real_data, axis=1))
Out[233]:
[<matplotlib.lines.Line2D at 0x7fb45dd18b90>]

Still not even numbers of bits per-packet.

In [238]: 
hist(np.sum(real_data, axis=1) % 2); None

And we don't seem to conserve parity in here. Let's see if that's generally true.

In [262]: 
bits_set = []

candidate_offsets = range(10, real_data.shape[1])

for i in candidate_offsets:
    parities = np.sum(real_data[:,0:i], axis=1) % 2
    n = np.sum(parities)
    bits_set.append(n)
    
plot(candidate_offsets, np.array(bits_set)/float(real_data.shape[0]))
ylim(0,1)
Out[262]:
(0, 1)

Well, doesn't look like any of those are simple parity bits, so this isn't a trivially-checksum'd packet.

What do the deltas look like here?

In [263]: 
imshow( (real_data[:-1,:] ^ real_data[1:,:]).T, interpolation="none", cmap=cm.gray_r)
Out[263]:
<matplotlib.image.AxesImage at 0x7fb45da60fd0>

Hrm. Let's reverse this.

In [229]: 
real_data_rev = real_data.copy()[:,::-1]
x = real_data_rev.tolist()
x.sort()
real_data_rev = np.array(x)
imshow(real_data_rev.T, interpolation="none", cmap=cm.gray_r)
Out[229]:
<matplotlib.image.AxesImage at 0x7fb45dfd3a10>

Sigh. I don't want to brute-force edit distance sort, but, here we go!

In [270]: 
distances = np.zeros((real_data.shape[0],real_data.shape[0]))
for i in range(real_data.shape[0]):
    for j in range(i+1, real_data.shape[0]):
        d = np.sum(real_data[i,:] ^ real_data[j,:])
        distances[i,j] = d
        distances[j,i] = d
        
hist(distances.flatten(), bins=range(real_data.shape[0])); None

See that empty 1 box there. That's proof that god hates us.

In [287]: 
to_show = []
to_show_h = []


for p in packets_peakish_slow:
    preamble_header = find_actual_preamble(np.array(p['bitstream']))
    if -1 == preamble_header:
        continue # preamble_header = 0
    to_add = p['bitstream'][(preamble_header):] + [0] * (125 - len(p['bitstream']) + preamble_header)
    to_show.append(to_add)

    to_show_h.append( hits[p['hit_i']] )
    

to_show = np.array(to_show)
real_data = to_show[:, 32+3:32+25]
In [288]: 
imshow(real_data.T, interpolation="none", cmap=cm.gray_r)
Out[288]:
<matplotlib.image.AxesImage at 0x7fb42ce55c10>
In [301]: 
idx = range(len(to_show_h))
idx.sort(key=lambda i: to_show_h[i]["Fc"])
In [302]: 
imshow(real_data[idx,:].T, interpolation="none", cmap=cm.gray_r)
Out[302]:
<matplotlib.image.AxesImage at 0x7fb45d179850>
In [309]: 
keys = ["Fc", "s_low", "s_length", "bw" ]

for i, k in enumerate(keys):
    idx.sort(key = lambda i: to_show_h[i][k])
    subplot(2,2,i+1)
    imshow(real_data[idx,:].T, interpolation="none", cmap=cm.gray_r)
    title("Sorted by " + k)
    
savefig("ordering_wtf.png", bbox="tight")

Okay, it looks like this is just somebody rolling through codes for some reason. Not sure what it is, but it's what it is.

Trying more 9600 baud material

In [314]: 
packets_slow = [ p for p in packets if p['baudrate'] > 9000 and p['baudrate'] < 10000 and len(p['bitstream']) > 150 ]

len(packets_slow)
Out[314]:
0

Whomp whomp. Nothing. On to the faster/bigger stuff.

In [320]: 
packets_fast = [ p for p in packets if p['baudrate'] > 50000 and p['baudrate'] < 200000 ]
hist([p['baudrate'] for p in packets_fast ], bins=100); None
In [409]: 
packets_fast = [ p for p in packets if p['baudrate'] > 90000 and p['baudrate'] < 110000 and len(p['bitstream']) > 200 ]
hist([p['baudrate'] for p in packets_fast ], bins=100); None
In [410]: 
to_show = []
to_show_h = []

l = 200

for p in packets_fast:
    to_add = p['bitstream'][:l] + [0] * (l - len(p['bitstream']))
    to_show.append(to_add)
    to_show_h.append( hits[p['hit_i']] )

    
to_show = np.array(to_show)
to_show.shape
Out[410]:
(331, 200)
In [411]: 
imshow(to_show, interpolation="none", cmap=cm.gray_r)
Out[411]:
<matplotlib.image.AxesImage at 0x7fb457b76b90>
In [412]: 
keys = ["Fc", "s_low", "s_length", "bw" ]
idx = range(len(to_show))
for i, k in enumerate(keys):
    idx.sort(key = lambda i: to_show_h[i][k])
    subplot(2,2,i+1)
    imshow(to_show[idx,:], interpolation="none", cmap=cm.gray_r)
    title("Sorted by " + k)
In [413]: 
to_show = []
to_show_h = []

l = 200

for p in packets_fast:
    bs = np.array(p['bitstream'])
    if len(bs) < 10:
        continue
    s = find_preamble_perfectly(bs)
    if None == s or s < 3 or s > l-10:
        continue
    to_add_head = list(p['bitstream'][s:s+l])
    to_add = to_add_head + [0] * (l - len(to_add_head))
    to_show.append(to_add)
    to_show_h.append( hits[p['hit_i']] )

to_show = np.array(to_show)
to_show.shape
Out[413]:
(327, 200)
In [419]: 
imshow(to_show.T, interpolation="none", cmap=cm.gray_r)
Out[419]:
<matplotlib.image.AxesImage at 0x7fb4579fd050>
In [415]: 
np.mean([ p['baudrate'] for p in packets_fast])
Out[415]:
99749.26008431938
In [420]: 
idx = range(len(to_show))
idx.sort(key = lambda i: to_show_h[i]['Fc'])
imshow(to_show[idx,:].T, interpolation="none", cmap=cm.gray_r)
title("Ordered by frequency")
Out[420]:
<matplotlib.text.Text at 0x7fb45796bfd0>
In [417]: 
idx = range(len(to_show))
idx.sort(key = lambda i: to_show_h[i]['s_low'])
imshow(to_show[idx,:], interpolation="none", cmap=cm.gray_r)
title("Ordered by transmission start time")
Out[417]:
<matplotlib.text.Text at 0x7fb457ba2c50>
In [418]: 
idx = range(len(to_show))
idx.sort(key = lambda i: to_show_h[i]['s_length'])
imshow(to_show[idx,:], interpolation="none", cmap=cm.gray_r)
title("Ordered by bitstream_length")
Out[418]:
<matplotlib.text.Text at 0x7fb457c0e850>
In [418]:
In [408]:
In [408]:
In []: